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' Abstract 

O 

We studied through Monte Carlo simulations, the kinetics of the two-species diffusion-limited 
reaction model with same species excluded volume interaction in substrates embedded on a 

in 

We study the time evolution of the concentration of single-particle species for various values 

o 



of substrate occupancies, 0.5927460, 0.61, 0.63, 0.65, 0.7, 0.8, and 1, where the first value 
corresponds to the percolating probability of the square lattice. We show that in the diffusion- 
limited reaction regime, the kinetics strongly depends on the presence of a bias along a particular 
square lattice direction, representing the net effect of a driving field. We were able to explain 

X 

the slow dynamics at high values of the driving field in terms of traps appearing in diluted 

a ; 

substrates, particularly at the percolation threshold. 
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I. INTRODUCTION 



Non-equilibrium phenomena are ubiquitous in nature and of great interest to model 



3, 



5|. Such a model of non-equili 



limited reaction model 



)del oi non-equili 
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jrium phenomena is the so-called diffusion- 
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12| . This model applies, when the time 



for two particles to react is less than the time for the particles to diffuse between colli- 
sion. In this limit, diffusion becomes the slowest, rate determining process. There has 
been a host of results in the last two decades regarding the various properties of these 
reactions, and the realm of application exceeds what the name suggests. In fact, diffusion- 
limited reactions, apart from an interest on their own as case study of non-equilibrium 
phenomena, have been proved useful in the understanding of a wide range of natural 
phenomena like colloids and aerosols dynamics, soliton-antisoliton dynamics in quasi- Id 
systems, vapor-de pos ited thin films, exciton reactions, and electron-hole recombination 



in semiconductors 
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The underlying reason for theoretical inter- 
est is related to the failure of the simple minded classical rate equation approach to the 
study of diffusion-limited phenomena. More specifically, the kinetics of these reactions 
is dominated by fluctuations in particle density, which cannot be grasped by a mean- 
field, rate-equation type of approach. For example, the bimolecular reaction A + B — > 
has an upper critical spatial dimension d c = 4, which means that below four spatial 
dimensions, the kinetics is dominated by fluctuations in particle density. As fluctua- 
tions are more effective on lower spatial dimensions, there are a host of results in d — 1 

, but some research has also been done in two and three dimen- 
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cases of the effect of biased diffusion 
field, on the kinetics o 
percolating cluster |fl, 



231 1 . In the present work, we 
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'ocus our attention to the relevant 



271 ]. as a consequence of a driving 



diffusion- limited reactions on substrates ranging from the diluted, 
to the compact substrate. We perform detailed analy- 
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sis of the various stages of evolution of the kinetics in the isotropic and biased diffusion 



cases, and for different substrate occupancies. We show novel interesting properties of 
these models including trap-escape dominated kinetics. 

The paper is organized as follows: In Section m we define the model and present the 
details of our simulations. Next, in Section IIII[ we present our preliminary results and 
discuss them. Finally, in Section HVj we present our conclusions. 



II. MODEL AND SIMULATIONS 



To simulate the model, we start by defining it in a square lattice with periodic boundary 
conditions for concreteness, but generalization to other lattices is straightforward. Accord- 
ing to the percolation model each lattice site is occupied with probability p and left empty 



with probability 1—p 



31 



32| . For values of the occupation probability p > p c , where p c is 
the percolation threshold probability, the largest cluster thus formed by nearest-neighbor 
sites has, at least, one path connecting opposite borders. In this percolative cluster, we 
deposit a uniform initial concentration of A- and B- type of reactive species, Cq, with equal 
number of both particle species. 

The kinetics of the system takes place on the percolating cluster defined above. Evo- 
lution takes place (Fig. H]a)) by randomly choosing a particle and allowing it to attempt 
to move to one of its nearest-neighbor sites. At every hopping attempt, only four events 
can take place. For example, a particle can diffuse by hopping to a nearest neighbor 
unoccupied site belonging to the cluster. As a result, in a diffusive move a particle solely 
changes its position (Fig. H]a)-i)). However, if the neighboring site is occupied by a like- 
species particle, excluded volume interaction prevents the particle from moving and the 
attempt fails, leaving both particles in their original positions (Fig. ^a)-ii)). A failed 
hopping attempt also takes place when the chosen nearest neighboring site belongs to the 
perimeter of the percolating cluster (Fig. [T] a)-iii)), and, as a result, the particle just re- 
mains in its original position. This is equivalent to having reflective boundary conditions 
within the percolating cluster. The last event occurs when the selected site hops onto an 
unlike-species particle. In this case, both particles react and desorb from the substrate 
(Fig. n a )-i v )) ; i-e., A + B — > 0. Finally, time is incremented, by the reciprocal of the 
remaining number of particles, whether an attempt is successful. 

In this work we simulated systems, all embedded in a square lattice with a size of 
1024 2 sites, with a substrate occupancy ranging from the percolation threshold given by 



p c = 0.5927460 up to the compact case with p — 1. The intermediate values of the 
substrate occupancy, p, we simulated are 0.61, 0.63, 0.65, 0.70, and 0.80. Applying an 
external field along one of the square lattice directions, causes the number of attempts 
of a particle to increase along this direction as compared to attempts in the opposite 
direction. The direction of the field is kept fixed along a square lattice direction, so, 
every n attempt in the field direction, a particle attempts to move, on average, pfii times 
along the field direction and p b n times in the opposite direction, where pf and pb are the 
probabilities of a particle to move, respectively, forward and backwards along the field 
direction (pf+Pb = 1)- We further define r\ = Pf—p b , which describes the bias introduced 
by the presence of the field. As defined, the bias r] — 0, corresponds to zero field condition, 
i.e., isotropic diffusion, while for rj = 1 we have a full bias condition. Perpendicularly to 
the field direction, diffusion is isotropic. Therefore, the probability, Vi, for a particle to 
attempt to move along any of the four nearest-neighbor directions of the square lattice is 
given by 

i±2 a i on g 

^ opposite to the field direction. (1) 

\ perpendicular to 

In the present work, the values of rj used in simulations are 0, 0.4, 0.8, and 0.9. 

The direct consequence of the presence of bias, is to favor the rapid build up of cor- 
relations along the direction of the field as compared, for example, to a perpendicular 
direction. Consequently, it is appropriate to use rectangularly shaped lattices, which take 
into account the above fact with the longer side of the substrate along the field, as ex- 
pected: we take, L\\/L± = (1 + rj), where Ly and Lj_ correspond to the lattice linear 
dimension along the field and perpendicular to it, respectively. 

A set of five samples were simulated for each pair of parameters, (p,rj), with an initial 
total concentration, c = 0.4, with periodic boundary conditions. 



III. RESULTS AND DISCUSSION 



We present our results for the time evolution of single-species concentration in two 
limiting cases, namely zero, rj = and strong ,r] = 0.95, driving fields for all the p values 
presented above. We start by presenting the time dependence of the concentration, as 
shown in Fig.|21 in the isotropic, zero field case. We observe three major distinct evolution 
stages (Fig. El- inset). The first stage lasts until the concentration drops by 18% at p c 
and by 27% in the compact system, and kinetics is well-characterized by rate equations. 
This is so, since the distribution of both particle species is homogeneous, thus spatial 
fluctuations are unimportant. Support for this line of reasoning comes from the fact that 
curves from different substrate occupancies differ only slightly within this stage, a clear 
indication of the lack of sensitivity on the details of the substrate. Moreover, for every 
simulation we also adopted a different substrate, so evolution is also insensitive to the 
substrate structure. Being on the first state decaying, the kinetics is now dominated by 
fluctuations on particle density, i.e., the diffusion-limited regime. As particles are depleted, 
i.e., annihilated during the first stage, domains start to form due to fluctuations in particle 
number. Regions rich in A-particles form a domain of such particles, while domains of 
.B-particles appear in regions with deficiency of A-particles. Moreover, particles have to 
diffuse until they collide, which represents the slowest, rate determining, process present 
in the system. The kinetics of these systems is profoundly related to the initial conditions, 
both spatial and temporal correlations grow over time due to local, initial fluctuations in 
particle- number of each species. The concentration decays as a power-law, t~ a , where the 
exponent a takes a value of 0.5, for the compact system in 2-D. In this regime, mean-field, 
classical rate equations do not describe the kinetics of the system well, since from these 
a concentration decay with a = 1 is obtained. Rate equations are not able to capture 
the spatial fluctuations in particle number, but these are fundamental for the proper 
characterization of the kinetics of diffusion-limited reaction systems. At the percolation 



threshold a ~ 0.347, which agrees with the result quoted in reference 3JJ, |3J|. The late 
times stage corresponds to a concentration lower than 10~ 3 , it is a consequence of the 
finiteness of the system, and it decays exponentially 

The various stages can be identified, including the value of a, by plotting the time 
dependence of the derivative of the concentration natural logarithm in order of the time 



natural logarithm shown in Figure |3] 34] . We plot the derivative versus time for all p 
values for the zero bias, isotropic diffusion case shown in Fig. |3J For t < 1, the various 
derivatives are not sensitive to the value of p, only spreading slightly at the end of the first 
stage. For t > 1, the system undergoes a transient stage until t ~ 10 3 time units before 
it reaches the diffusion-limited regime for all values of p > 0.65, while for p = p c , 0.61 the 
transient time is ~ 50. The differences in the kinetics for different values of p stems from 
the fact that the substrates geometrical properties change from fractal, scale-invariant 



behavior at p c to a translational invariant 



2 1, and also from the domain formation due 



to local density fluctuations in particle number. Time spent in the transient stages may 
last long enough to prevent a significantly extended diffusion-limited regime from being 
observed, even for a lattice of size 1024 2 sites, before late stages kinetics kicks in, which 
shows difficulties simulating these systems. 

Regarding the strong bias case, the differences relative to the isotropic case emerge right 
after the first stage as shown in Fig. 0] for all values of p, but the compact system, where 
the overall behavior follows the unbiased diffusion case. In contrast to the isotropic case, 
the transition from the transient to the diffusion-limited regime is now more incisive, 
and the transition times also vary with p. Therefore, the effect of a strong bias is to 
enhance the effect of topological restrictions of the substrate on the kinetics. Such an 
effect of the bias is to slow down the kinetics, a counter intuitive behavior, since the 
presence of a bias generates a net drift of particles along its direction. As can be seen 
from Fig El the curve for p — 0.8 shows faster decay than the isotropic case over the 



transient regime. These differences are confirmed by the presence of a well defined peak 
in the transient stage, as shown in Fig. 03 followed by flat terraces. The height and width 
of the peak of time derivative, see Fig. now also depends on substrate occupancy. 
For substrate occupancies, values of p ~ 1, substrate imperfections are comparable to 
obstacles. However, for occupancies values of p near p c , the substrate is constituted by 
topological traps with few or no exits in the direction of the bias, so particles take a long 
time to escape from such traps. Moreover, in the diffusion-limited regime, the exponent 
rapidly converges to a value, different from its zero bias counterpart for all values of p, 
but p = 1 [35 1 . 

We finalize the presentation of our results by presenting the concentration time evolu- 
tion for various values of the applied bias for a substrate occupancy of 0.65, as shown in 
Figure [HI We observe with increasing values of r], i.e., of the bias, a clear slowing down 
of the kinetics. The slow down of the kinetics is so accentuated that simulations cannot 
last long enough to observe late stage kinetics. In Fig. we show the time dependence of 
the numerical derivative, as defined above, for various values of rj, namely, 0, 0.4, 0.8, and 
0.9. Notice, the wide variation of late stage kinetic exponents, keeping in mind that for 
values of rj > 0, one can only make rough estimates, since we did not enter deep enough 
into the late stage kinetics. Finally, we present in Table HJ a summary of the values of the 
exponent a. 



IV. CONCLUSION 



Though non-compact substrates are quite common in nature, little effort has been put 
on the characterization of effects of topological restrictions on the kinetics. We performed 
a series of simulations on a square lattice with 1024 2 sites for various values of substrate 
occupancy ranging from the percolation threshold to a compact system. Our preliminary 
results show non-trivial substrate occupancy effects on the various stages of the kinetics. 
Further inclusion of bias, i.e., preferred diffusion along one the directions of a driving field, 
enhances the topological constraints present in the substrate. One observes substantial 
slow down of the kinetics for values of the bias, 77, above 0.4. Substrate topological 
restrictions, together with the presence of a bias, enhance early times reactivity, but 
slows it down at long times. Our results can be of relevance to researchers working on 
heterogenous reactions on surfaces and more generally to researchers working on non- 
compact substrates problems. 

Acknowledgments 

This research has been partially supported by PRODEP III (Application 
1/5.2/PRODEP/1996 University Advanced Education - Student 10 (NJG)) and Fundacao 
para a Ciencia e a Tecnologia (grant POCTI/CTM/41574/2001 - Computational Nano- 
physics) . 



[1] J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge 

University Press, Cambridge, 1999). 
[2] A. Bunde and S. Havlin, eds., Fractals and Disordered Systems (Springer, Berlin, 1991). 
[3] D. ben-Avraham and S. Havlin, Diffusion and Reactions in Fractals and Disordered Systems 

(Cambridge University Press, Cambridge, 2000). 



[4] V. Privman, Trends in Stat. Phys. 1, 89 (1994). 

[5] V. Privman, ed., Nonequilibrium Statistical Mechanics in One Dimension (Cambridge Uni- 
versity Press, Cambridge, 1997). 

[6] A. A. Ovchinnikov and Y. B. Zeld'ovitch, Chem. Phys. 28, 215 (1972). 

[7] D. C. Torney and H. M. McConnell, Proc. R. Soc. London, Ser. A 387, 147 (1983). 

[8] D. Toussaint and F. Wilczek, J. Chem. Phys. 78, 2642 (1983). 

[9] K. Kang and S. Redner, Phys. Rev. Lett. 52, 955 (1984). 

[10] K. Kang and S. Redner, Phys. Rev. A 30, 2833 (1984). 

[11] V. Kuzovkov and E. Kotomin, Rep. Prog. Phys. 51, 1479 (1988). 

[12] A. Cadilhe, M. L. Glasser, and V. Privman, Int. J. Mod. Phys. B 11, 109 (1997). 

[13] D. ben-Avraham, M. A. Burschka, and C. R. Doering, J. Stat. Phys. 60, 695 (1990). 

[14] S. K. Friedlander, ed., Smoke Dust and Haze: Fundamentals of Aerosol Behavior (John 

Wiley, New York, 1977). 

[15] F. Family and P. Meakin, Phys. Rev. A 40, 3836 (1989). 

[16] R. Kopelman, Phil. Mag. B 56, 717 (1987). 

[17] R. Kopelman, S. J. Parus, and L. Prasad, Chem. Phys. 128, 209 (1988). 

[18] E. A. Schiff, J. Non-Cryst. Solids 190, 1 (1995). 

[19] D. C. Torney and H. M. McConnell, J. Phys. Chem. 87, 1941 (1983). 

[20] Z. Racz, Phys. Rev. Lett 55, 1707 (1985). 

[21] M. Bramson and J. L. Lebowitz, Phys. Rev. Lett. 61, 2397 (1988). 

[22] V. Privman, J. Stat. Phys. 69, 629 (1992). 

[23] D. ben-Avraham and C. R. Doering, Phys. Rev. A 37, 5007 (1988). 

[24] S. A. Janowsky, Phys. Rev. E 52, 2535 (1995). 

[25] S. A. Janowsky, Phys. Rev. E 51, 1858 (1995). 

[26] I. Ispolatov, P. Krapivsky, and S. Redner, Phys. Rev. E 52, 2540 (1995). 



[27] Y. Shafrir and D. ben-Avraham, Phys. Lett. A 278, 184 (2001). 
[28] P. Meakin and H. E. Stanley, J. Phys. A: Gen. Math. 17, L173 (1984). 
[29] L. W. Anacker and R. Kopelman, J. Phys. Chem. 91, 5555 (1987). 
[30] P. Argyrakis and R. Kopelman, Phys. Rev. A 45, 5814 (1992). 
[31] D. Stauffer, Phys. Rep. 54, 1 (1979). 

[32] D. Stauffer and A. Aharony, Introduction to Percolation Theory (Taylor & Francis, London, 
1994), 2nd ed. 

[33] N. J. Goncalves, J. M. A. S. Duarte, and A. Cadilhe, Preprint (2004). 
[34] N. J. Gongalves, J. M. A. S. Duarte, and A. Cadilhe, Preprint (2004). 
[35] More extended accounts of this line of reasoning and of the exponent values in the ther- 
modynamic limit will be published in a different article, due to limitations in article size. 



TABLE I: Exponent a in the diffusion-limited regime for each simulated pair of substrate occu- 
pancy, p, and anisotropy, rj. 



V 

p 0.0 0.4 0.8 0.9 

p c 0.36 ± 0.06 0.07 ± 0.02 0.06 ± 0.01 0.05 ± 0.01 

0.61 0.48 ± 0.04 0.10 ± 0.01 0.06 ± 0.01 0.06 ± 0.01 

0.63 0.51 ± 0.05 0.16 ± 0.02 0.08 ± 0.01 0.07 ± 0.01 

0.65 0.55 ± 0.05 0.24 ± 0.03 0.11 ± 0.01 0.09 ± 0.01 

0.70 0.52 ± 0.04 0.58 ± 0.04 0.17 ± 0.04 0.17 ± 0.02 

0.80 0.55 ± 0.04 0.64 ± 0.07 0.50 ± 0.17 0.27 ± 0.10 

1 0.52 ± 0.07 0.57 ± 0.04 0.58 ± 0.04 0.54 ± 0.03 



FIG. 1: Rules of the A + B— > diffusion-limited reaction with same species excluded volume 
interaction on fractal percolating clusters, a) Typical configuration at early times. Dark grey 
regions represent the percolating cluster, light grey regions represent remaining clusters, and 
white regions represent empty sites. The image represents a cutout of a larger system of 1024 2 
sites, b) From the local configuration, highlighted in a), there is a i?-particle for which a hopping 
attempt can lead to: i) a diffusion - upwards; ii) an excluded volume - rightwards; iii) bouncing 
off from the boundary - downwards; iv) a reaction - leftwards. In all cases i)-iv) light grey 
squares represent next-nearest neighbors, which do not directly contribute to the kinetics. 

FIG. 2: Single-species concentration versus time log-log plot at zero bias for various substrate 
occupancies, p, on a 1024 2 sites system. 

FIG. 3: Representation of the local slopes of — d(logc(t))/d(logt), with r] = 0, for various 
substrate occupancies. 

FIG. 4: Single-species concentration versus time log-log plot at a bias of 0.9 for various substrate 
occupancies, p, on a 1024 2 sites system. 

FIG. 5: Representation of the local slopes of — d(logc(t))/d(logt), with r] = 0.9, for various 
substrate occupancies. 

FIG. 6: Single-species concentration versus time log-log plot at zero bias for a substrate occu- 
pancies of 0.65, on a 1024 2 sites system. 

FIG. 7: Representation of the local slopes of — d(logc(t))/d(logt), for various values of 77, and 
for a substrate occupancy of 0.65. 
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